module vec_mod; 
use dbg_mod; 
implicit none;
!private;
logical, save :: DBG_VEC_NEW = .TRUE., SAVE_VEC_NEW;
logical, save :: DBG_VEC_IF = .TRUE., SAVE_VEC_IF;
logical, save :: DBG_VEC_ENTRY = .FALSE., SAVE_VEC_ENTRY;
!logical, save :: DBG_VEC_ENTRY = .TRUE.;
integer, parameter :: DUFF_SIZE = 8;
integer, parameter :: DUFF_SIZE_16 = 16;
integer, save :: NEW_VEC_COUNT = 0;

contains; 

subroutine dbg_vec_to_save();  implicit none;  
  SAVE_VEC_NEW = DBG_VEC_NEW;
  SAVE_VEC_IF = DBG_VEC_IF;
  SAVE_VEC_ENTRY = DBG_VEC_ENTRY;
end subroutine
subroutine dbg_vec_from_save();  implicit none;  
  DBG_VEC_NEW = SAVE_VEC_NEW;
  DBG_VEC_IF = SAVE_VEC_IF;
  DBG_VEC_ENTRY = SAVE_VEC_ENTRY;
end subroutine
subroutine dbg_vec_all_on();  implicit none;  
  DBG_VEC_NEW = .TRUE.;
  DBG_VEC_IF = .TRUE.;
  DBG_VEC_ENTRY = .TRUE.;
end subroutine

function new_vec(n) result(vec_p);    implicit none;  
  integer, intent(in) :: n;  integer :: status;
  real(R14), dimension(:), pointer :: vec_p;      
  allocate(vec_p(n), STAT=status);
  call dbg_info("-->new_vec()", DBG_VEC_ENTRY);
  vec_p = 0.;     call dbg_vec("ptr=", vec_p);
  NEW_VEC_COUNT = NEW_VEC_COUNT + 1;
  call dbg_int("NEW_VEC_COUNT=", NEW_VEC_COUNT, DBG_VEC_NEW);
  return;
end function

function new_vec_from(from_p) result(vec_p);   implicit none;  
  integer :: status, min_idx, max_idx;
  real(R14), dimension(:), pointer :: from_p, vec_p;      
  call dbg_info("-->new_vec_from()", DBG_VEC_ENTRY);
  min_idx = LBOUND(from_p, 1);
  max_idx = UBOUND(from_p, 1);
  allocate(vec_p(min_idx:max_idx), STAT=status);
  vec_p = 0.;       call dbg_vec("ptr=", vec_p);
  NEW_VEC_COUNT = NEW_VEC_COUNT + 1;
  call dbg_int("NEW_VEC_COUNT=", NEW_VEC_COUNT, DBG_VEC_NEW);
  return;
end function

subroutine del_vec(pv);  implicit none;  
  real(R14), dimension(:), pointer :: pv;      
  call dbg_info("-->del_vec()", DBG_VEC_ENTRY);
  call dbg_bool("ASSOCIATED(pv)=", ASSOCIATED(pv), DBG_VEC_IF);
  if (ASSOCIATED(pv)) then
    DEALLOCATE(pv);    !NULLIFY(ptr); !DEALLOCATE() makes it null
    NEW_VEC_COUNT = NEW_VEC_COUNT - 1;
    call dbg_int("NEW_VEC_COUNT=", NEW_VEC_COUNT, DBG_VEC_NEW);
  end if     
  return;
end subroutine

subroutine del_vec2(pv, pv2);  implicit none;  
  real(R14), dimension(:), pointer :: pv, pv2;      
  call del_vec(pv);
  call del_vec(pv2);
  return;
end subroutine

subroutine del_vec3(pv, pv2, pv3);  implicit none;  
  real(R14), dimension(:), pointer :: pv, pv2, pv3;      
  call del_vec(pv);
  call del_vec(pv2);
  call del_vec(pv3);
  return;
end subroutine

real(R14) function vec_dot_prod(v, v2) result(res);  implicit none;  
  real(R14), dimension(:), intent(in) :: v, v2;  
  call dbg_info("-->vec_dot_prod()", DBG_VEC_ENTRY);
  res = DOT_PRODUCT(v, v2);
  return
end function

!real(R14) function vec_dot_prod2(v, v2) result(res);  implicit none;  
!  real(R14), dimension(:), intent(in) :: v, v2;  
!  call dbg_info("-->vec_dot_prod2()", DBG_VEC_ENTRY);
!  res = SUM(v * v2);
!  return
!end function

real(R14) function vec_dot_prod3(v, v2, v3) result(res);  implicit none;  
  real(R14), dimension(:), intent(in) :: v, v2, v3;  
  call dbg_info("-->vec_dot_prod3()", DBG_VEC_ENTRY);
  res = SUM(v * v2 * v3);
  return
end function

real(R14) function vec_dot_prod4(v, v2, v3, v4) result(res);  implicit none;  
  real(R14), dimension(:), intent(in) :: v, v2, v3, v4;  
  call dbg_info("-->vec_dot_prod4()", DBG_VEC_ENTRY);
  res = SUM(v * v2 * v3 * v4);
  return
end function

!================================ DO NOT USE BELOW ============

real(R14) function vec_dot_do(v, v2) result(res);  implicit none;  
  real(R14), dimension(:), intent(in) :: v, v2;   
  integer :: i;
  if (DBG_VEC_ENTRY) write(*,*) "vec_dot_do";
  res = 0.;
  do i = LBOUND(v, 1), UBOUND(v, 1)
    res = res + v(i) * v2(i);
  end do 
  return
end function

real(R14) function vec_dot_do3(v, v2, v3) result(res);  implicit none;  
  real(R14), dimension(:), intent(in) :: v, v2, v3;   
  integer :: i;
  if (DBG_VEC_ENTRY) write(*,*) "vec_dot_do";
  res = 0.;
  do i = LBOUND(v, 1), UBOUND(v, 1)
    res = res + v(i) * v2(i) * v3(i);
  end do 
  return
end function


real(R14) function vec_dot_duff(v, v2) result(res);  implicit none;  
  real(R14), dimension(:), intent(in) :: v, v2;   
  integer :: i, m, i8, min_i, max_i, n, start_m, start_i8;
  if (DBG_VEC_ENTRY) write(*,*) "vec_dot_duff";
  min_i = LBOUND(v, 1); max_i = UBOUND(v, 1);
  res = 0.;
  i = min_i;
  n = max_i - min_i + 1;         !call dbg_int("n=", n);
  start_m = MOD(n, DUFF_SIZE);   !call dbg_int("start_m=", start_m);
  do m = start_m, 1, -1 
    res = res + v(i) * v2(i);  i = i + 1;      
    !call dbg_int("i=", i);
  end do
  start_i8 = n / DUFF_SIZE;
  !call dbg_int("start_i8=", start_i8);
  !call dbg_int("before do i=", i);
  do i8 = start_i8, 1, -1
    res = res + v(i) * v2(i);      i = i + 1; !1
    res = res + v(i) * v2(i);      i = i + 1; !2
    res = res + v(i) * v2(i);      i = i + 1; !3
    res = res + v(i) * v2(i);      i = i + 1; !4
    res = res + v(i) * v2(i);      i = i + 1; !5
    res = res + v(i) * v2(i);      i = i + 1; !6
    res = res + v(i) * v2(i);      i = i + 1; !7
    res = res + v(i) * v2(i);      i = i + 1; !8
  end do
  !call dbg_int("after do i=", i);
  return
end function

real(R14) function vec_dot_duff3(v, v2, v3) result(res);  implicit none;  
  real(R14), dimension(:), intent(in) :: v, v2, v3;   
  integer :: i, m, i8, min_i, max_i, n, start_m, start_i8;
  if (DBG_VEC_ENTRY) write(*,*) "vec_dot_duff";
  min_i = LBOUND(v, 1); max_i = UBOUND(v, 1);
  res = 0.;
  i = min_i;
  n = max_i - min_i + 1;         !call dbg_int("n=", n);
  start_m = MOD(n, DUFF_SIZE);   !call dbg_int("start_m=", start_m);
  do m = start_m, 1, -1 
    res = res + v(i) * v2(i)* v3(i);  i = i + 1;      
    !call dbg_int("i=", i);
  end do
  start_i8 = n / DUFF_SIZE;
  !call dbg_int("start_i8=", start_i8);
  !call dbg_int("before do i=", i);
  do i8 = start_i8, 1, -1
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !1
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !2
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !3
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !4
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !5
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !6
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !7
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !8
  end do
  !call dbg_int("after do i=", i);
  return
end function

real(R14) function vec_dot_duff3c(v, v2, v3) result(res);  implicit none;  
  real(R14), dimension(:), intent(in) :: v, v2, v3;   
  integer :: i, m, i8, min_i, max_i, n, start_m, start_i8;
  if (DBG_VEC_ENTRY) write(*,*) "vec_dot_duff";
  min_i = LBOUND(v, 1); max_i = UBOUND(v, 1);
  res = 0.;
  i = min_i;
  n = max_i - min_i + 1;         !call dbg_int("n=", n);
  start_m = MOD(n, DUFF_SIZE_16);   !call dbg_int("start_m=", start_m);
  do m = start_m, 1, -1 
    res = res + v(i) * v2(i)* v3(i);  i = i + 1;      
    !call dbg_int("i=", i);
  end do
  start_i8 = n / DUFF_SIZE_16;
  !call dbg_int("start_i8=", start_i8);
  !call dbg_int("before do i=", i);
  do i8 = start_i8, 1, -1
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !1
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !2
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !3
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !4
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !5
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !6
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !7
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !8
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !9
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !10
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !11
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !12
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !13
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !14
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !15
    res = res + v(i) * v2(i) * v3(i);      i = i + 1; !16
  end do
  !call dbg_int("after do i=", i);
  return
end function

real(R14) function vec_dot_duff3b(v, v2, v3) result(res);  implicit none;  
  real(R14), dimension(:), intent(in) :: v, v2, v3;   
  integer :: i, m, i8, min_i, max_i, n, start_m, start_i8;
  if (DBG_VEC_ENTRY) write(*,*) "vec_dot_duff";
  min_i = LBOUND(v, 1); max_i = UBOUND(v, 1);
  res = 0.;
  i = min_i;
  n = max_i - min_i + 1;         !call dbg_int("n=", n);
  start_m = MOD(n, DUFF_SIZE_16);   !call dbg_int("start_m=", start_m);
  do m = start_m, 1, -1 
    res = res + v(i) * v2(i)* v3(i);  i = i + 1;      
    !call dbg_int("i=", i);
  end do
  start_i8 = n / DUFF_SIZE_16;
  !call dbg_int("start_i8=", start_i8);
  !call dbg_int("before do i=", i);
  do i8 = start_i8, 1, -1
    res = res + v(i) * v2(i) * v3(i);      !i = i + 1; !1
    res = res + v(i+1) * v2(i+1) * v3(i+1);      !i = i + 1; !2
    res = res + v(i+2) * v2(i+2) * v3(i+2);      !i = i + 1; !3
    res = res + v(i+3) * v2(i+3) * v3(i+3);      !i = i + 1; !4
    res = res + v(i+4) * v2(i+4) * v3(i+4);      !i = i + 1; !5
    res = res + v(i+5) * v2(i+5) * v3(i+5);      !i = i + 1; !6
    res = res + v(i+6) * v2(i+6) * v3(i+6);      !i = i + 1; !7
    res = res + v(i+7) * v2(i+7) * v3(i+7);      
    res = res + v(i+8) * v2(i+8) * v3(i+8);      
    res = res + v(i+9) * v2(i+9) * v3(i+9);      
    res = res + v(i+10) * v2(i+10) * v3(i+10);      
    res = res + v(i+11) * v2(i+11) * v3(i+11);      
    res = res + v(i+12) * v2(i+12) * v3(i+12);      
    res = res + v(i+13) * v2(i+13) * v3(i+13);      
    res = res + v(i+14) * v2(i+14) * v3(i+14);      
    res = res + v(i+15) * v2(i+15) * v3(i+15);      
    i = i + DUFF_SIZE_16; 
  end do
  !call dbg_int("after do i=", i);
  return
end function
end module 